Certain types of simulators, which simulate processes such as chemical reactions, have operated by solving differential reaction-rate equations. See, for instance, the introductory portion of Gillespie, "Exact Stochastic Simulation of Coupled Chemical Reactions," Journal of Physical Chemistry, vol. 81, page 2340 (1977). Because such equations usually cannot be solved analytically, computerized numerical approximations have been used. However, such simulations are based on the not-always-realistic premises that the reactions are continuous and deterministic. First, computerized numerical solutions of such differential equations are inherently based on discrete points in time, rather than on a continuous time flow. Also, in the context of simulation through solution of such differential reaction rate equations, Gillespie describes "deterministic" as being an unrealistic condition, because molecular populations can only change by discrete integer amounts, and because the exact positions and velocities of the individual molecules, which must be taken into account for a fully deterministic analysis, are not available.
By contrast, stochastic simulators have been used to make discrete time simulations of chemical processes. Stochastic simulators operate by tracing a set of individual events which may occur over a time period during which the process is to be simulated. The simulators operate in a stepwise fashion, in terms of an independent state variable. The independent state variable is typically time, and the simulation operates over a given domain of values for the independent state variable, such as from an initial time to a conclusion time, or from an initial state at an initial time until a final state is reached.
The simulator makes a determination that each of the events has a certain probability of occurrence, which is determined by the system conditions and material properties, at any specific point of time. Typically, the probability is a function of the concentrations of the reactants, and of other factors such as temperature, pressure, etc. In a predefined time-window or interval, which is commonly referred to as a time step, any one of the simulated events may occur. There is, on the average, one occurrence of an event within a time step when the length of that time step is selected to be equivalent to the inverse of the total probability of occurrence of all events. The likelihood of occurrence of a certain event within that time step is proportional to the rate for that individual event.
A stochastic simulator commonly comprises at least one processor which is programmed to define a time step, based on the sum total of the reaction rates for all of the reactions being simulated, and then to randomly select one of those reactions, for which there will be an event which is to take place in that time step. The likelihood of an event being randomly selected is proportional to the probability of occurrence of the selected event, relative to that of the other events. As the consequence of the occurrence of the selected event, the processor updates the system conditions by use of material property tables, system process equations and other required data groups stored in at least one data storage apparatus. The system update may change the probabilities of occurrence of the events, resulting in a different total probability and a different time step. The processor then selects a new time step and repeats the event selection and system condition update processes, to propagate a system through the entire time span for system simulation.
A stochastic simulator was first used to simulate a spatially homogeneous chemical system where many chemical reactions occur (D. L. Bunker et al., Combustion and Flame, vol. 23, p. 373, 1974). The probability of an individual chemical reaction is dependent on the rate constants, and on the concentrations of a variety of chemical reactants. The functional dependence and the concentrations are stored in the storage apparatus. The processor generates a random number between 0 and 1 to select among various chemical reactions a specific reaction in each time step, thus allowing the simulation to propagate in time. The time base for determining the length of time step is based on the total sum of all reaction probabilities.
Additional pioneering work was done by Bunker and Houle in "User's Manual: MSIM4--Discrete Simulator for Kinetic Mechanisms," Univ. of California, Irvine, Jul. 12, 1974, and by Hinsberg and Houle in "MSIMPC v2.0: An Interactive Discrete Chemical Mechanism Simulator for the IBM PC (Laboratory and Classroom Version)," IBM RJ 7814 (Nov. 8, 1990), to implement a simulator similar to that of Bunker and Houle for the IBM Personal Computer.
Another chemical simulator has been developed by Gillespie wherein, for each stochastic time step, the simulator generates two random numbers, each between 0 and 1. The first random number is used for selecting a chemical reaction event which is to be simulated. The probabilities of all of the possible reactions are normalized from 0 to 1. Thus, the range from 0 to 1 is divided up into subintervals, each subinterval being assigned to one of the reactions, and having a size proportional to that reaction's probability. The first random number falls within the subinterval of one of the reactions. That reaction is selected as the event.
The second random number is used for properly weighting the time base used for time step determination. The reciprocal of the sum of the reaction rate, which is a time interval, is used to calculate a time step .tau. for the simulation, using a formula such as ##EQU1## where rate.sub.j is the reaction rate of the j-th reaction of the simulation, and R is the second random number between 0 and 1. It will be seen that, for such stochastic simulations, time steps need not be uniform in size, or predictable. (Gillespie, "A General Method for Numerically Simulating the Stochastic Time Evolution of Coupled Chemical Reactions," Journal of Computational Physics, Vol. 22, No. 4, pp. 403-434, 1976). Better accuracy is achieved with this method than with Bunker's method, when fluctuating systems are involved.
In both the Bunker and the Gillespie simulators, the probability for each possible reaction is determined from the rate of that reaction, and these rates are added to produce a total reaction rate, or, equivalently, a total probability of any reaction taking place. The reaction probabilities may be thought of as normalized reaction rates, since the probability has a value ranging from 0 to 1. A reaction probability P.sub.i for an i-th reaction in a system in which j number of reactions are to take place, may be given by the expression ##EQU2## where the rate.sub.i and rate.sub.j expressions are the non-normalized reaction rates.
The Gillespie and Bunker references used the terms "rate" and "probability" more or less interchangably. However, in the present specification, a rigorous distinction will be made between "rate" and "probability," to explain clearly how each relates to the stochastic simulations being described.
Conventional stochastic simulators have a drawback, however, related to the issue of computational efficiency. Frequently, a system to be simulated, which includes a set of chemical reactions, will include one or more reversible reactions. If a reaction is reversible, the simulation may encounter conditions under which the reversible reaction or reactions are in equilibrium. Moreover it will sometimes be the case that a reaction in equilibrium will also have a reaction probability, in either direction, much greater (such as, several orders of magnitude greater) than the other, non-equilibrium, reactions in the simulation. Thus, in long stretches of computer simulation time, the reactants shuffle back and forth in the two complementary equilibrium steps, and only rarely move the simulation toward completion by undergoing the non-equilibrium reaction. The processing time in between successive occurrences of the non-equilibrium reaction is, in effect, wasted, because it does not produce a material change in the state of the system being simulated. Simulations in which such equilibrium reactions occur have been of limited use in solving many problems of scientific and technological interest, because the calculations are simply too slow to be of practical use.
The Bunker and Houle MSIM4 simulator, referred to above, attempts to overcome this computational efficiency problem by using a pattern recognition routine which identifies steps which have entered equilibrium, halts the simulation, and sends an interactive message to the system operator. The user then circumvents the equilibrium by combining the equilibrium step with other steps of the simulated process, to avoid explicit simulation of the equilibrium step. One way to do this is to add together the stoichiometries of the equilibrium step and the non-equilibrium step, and then to multiply the rate constant by the equilibrium constant. However, this approach to solving the computational efficiency problem has certain limitations. First, it requires prior knowledge of the system to be simulated, and this knowledge may not be available to the user. Second, it cannot be applied in all possible simulations, but only to those simulations with appropriate preconditions. Third, systems in which transient equilibrium conditions occur only during certain time segments of the process require each segment of the process, with differing behavior, to be analyzed and solved independently, resulting in inefficient and inconvenient operation.